function plot_basis_function_3d(basis_type, flag)

N = reference_basis_function_3d(basis_type);

if nargin == 1
    Nb = size(N,1);
else
    Nb = min(size(N,1), flag);
end

switch basis_type
    case "P1"
        [x, y, z] = meshgrid(linspace(-1, 1, 100), linspace(-1, 1, 100), linspace(-1, 1, 100));
        for i = 1:Nb
            phi = N{i}(x, y, z);
            isosurface(x, y, z, phi);
            hold on;
        end
    
    otherwise
        error("Invalid basis type");
end


end